Attribution of Space‐Time Variability in Global‐Ocean Dissolved Inorganic Carbon

Abstract The inventory and variability of oceanic dissolved inorganic carbon (DIC) is driven by the interplay of physical, chemical, and biological processes. Quantifying the spatiotemporal variability of these drivers is crucial for a mechanistic understanding of the ocean carbon sink and its future trajectory. Here, we use the Estimating the Circulation and Climate of the Ocean‐Darwin ocean biogeochemistry state estimate to generate a global‐ocean, data‐constrained DIC budget and investigate how spatial and seasonal‐to‐interannual variability in three‐dimensional circulation, air‐sea CO2 flux, and biological processes have modulated the ocean sink for 1995–2018. Our results demonstrate substantial compensation between budget terms, resulting in distinct upper‐ocean carbon regimes. For example, boundary current regions have strong contributions from vertical diffusion while equatorial regions exhibit compensation between upwelling and biological processes. When integrated across the full ocean depth, the 24‐year DIC mass increase of 64 Pg C (2.7 Pg C year−1) primarily tracks the anthropogenic CO2 growth rate, with biological processes providing a small contribution of 2% (1.4 Pg C). In the upper 100 m, which stores roughly 13% (8.1 Pg C) of the global increase, we find that circulation provides the largest DIC gain (6.3 Pg C year−1) and biological processes are the largest loss (8.6 Pg C year−1). Interannual variability is dominated by vertical advection in equatorial regions, with the 1997–1998 El Niño‐Southern Oscillation causing the largest year‐to‐year change in upper‐ocean DIC (2.1 Pg C). Our results provide a novel, data‐constrained framework for an improved mechanistic understanding of natural and anthropogenic perturbations to the ocean sink.


Introduction
Since the outset of the industrial era, human activities such as fossil fuel combustion and cement production have rapidly increased the concentration of atmospheric carbon dioxide (CO 2 ). On timescales relevant to this anthropogenic perturbation, the ocean constitutes the largest non-geological carbon reservoir on Earth (37,200 ± 195 Pg C of DIC;Keppler et al., 2020), storing roughly 17 and 50 times more carbon than the terrestrial biosphere and pre-industrial atmosphere, respectively (Friedlingstein et al., 2019). As such, the ocean plays a critical role in regulating global climate (Archer et al., 2009) and mitigating anthropogenic warming by sequestering atmospheric CO 2 (i.e., the "ocean carbon sink"; Heinze et al., 2015). To date, the ocean sink has absorbed roughly 40% of industrial-era fossil carbon emissions (Khatiwala et al., 2009(Khatiwala et al., , 2013McKinley et al., 2017). This has motivated numerous studies to monitor patterns and trends in the ocean sink, with recent global estimates of anthropogenic carbon uptake yielding values of roughly 2.3 ± 0.6 (2000Khatiwala et al., 2009Khatiwala et al., ) and 2.6 ± 0.3 (1994Khatiwala et al., -2007Gruber et al., 2019) Pg C year −1 . A major outstanding question is whether the oceans will continue to absorb roughly 40% of anthropogenic CO 2 emissions in the future, or whether this fraction will change significantly as emissions and/or the climate state change.
The trajectory of the ocean sink is driven by processes that drive CO 2 undersaturation (ocean uptake) and oversaturation (ocean outgassing) in surface waters. Many studies suggest an increase in the ocean sink during recent decades (Sarmiento & Gruber, 2002;Khatiwala et al., 2013;DeVries, 2014;DeVries et al., , 2019; however, storage rates are subject to considerable variability over a range of space-time scales. This hinders our ability to predict the long-term variability and strength of the ocean sink (Landschützer et al., 2016;McKinley et al., 2020;Randerson et al., 2015) and limits efforts that inform policy in response to anthropogenic climate change (Dunne, 2016). Quantifying the magnitude and space-time variability of the ocean sink has been recognized as an important goal in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2013). In this study, we specifically focus on DIC as it: (a) is the sum of aqueous inorganic carbon species and represents ∼98% of the carbon in the ocean (Hansell, 2013) and (b) is directly linked to atmospheric CO 2 concentrations.
Trends and patterns in the ocean carbon reservoir, specifically the DIC pool, are driven by the complex interplay of physical, chemical, and biological processes. Diagnosing and quantifying the space-time variability of these drivers (Lauderdale et al., 2016;Lovenduski et al., 2013) is critical for a mechanistic understanding of the ocean DIC state and the future trajectory of its sources (DIC gain) and sinks (DIC loss). At the ocean surface, DIC concentration is modulated by ocean-atmosphere CO 2 exchange, which depends on CO 2 solubility, air-sea gas transfer velocity, and the difference between atmospheric and surface-ocean partial pressure of CO 2 (apCO 2 and pCO 2 , respectively). Additionally, natural and anthropogenic vertical DIC gradients in the water column exhibit opposing signs (i.e., natural DIC typically increases with depth while anthropogenic DIC decreases with depth).
Three-dimensional ocean circulation and mixing (DeVries et al., 2019;Silvano, 2020) subducts surface-ocean DIC to depth, which isolates it from the atmosphere, and also upwells deep remineralized DIC back to the surface. This ventilated DIC is further transformed by exchange with the atmosphere and upper and intermediate waters, eventually being subducted back to depth (Talley, 2013). Biological processes redistribute DIC by uptake and fixation, export of particulates, remineralization of dissolved and particulate organic matter, dissolution of particulate inorganic carbon, and sedimentation of particulate matter -collectively termed the "biological pump" (Boyd et al., 2019).
To better characterize the ocean sink response to natural and forced variability, studies have used statistical extrapolation of surface-ocean pCO 2 and gas transfer parameterizations (e.g., Wanninkhof, 1992) to reconstruct maps of air-sea CO 2 flux and its space-time variability (Landschützer et al., 2013;Rödenbeck et al., 2013Rödenbeck et al., , 2015Takahashi et al., 2002;Wanninkhof et al., 2013). Most of these extrapolations leverage statistical relationships between in-situ surface-ocean pCO 2 , remotely sensed sea-surface temperature (SST) and Chlorophyll (Chl), and climatological mixed-layer depth (MLD) and sea-surface salinity (SSS). Recent efforts have further extended these methods to monthly climatological DIC (Keppler et al., 2020) and surface-ocean carbonate chemistry parameters (Gregor & Gruber, 2021). However, these extrapolations suffer from a number of methodological uncertainties, data sparsity, and observation noise . For example, since the mid-1980s, less than 2% of all monthly 1° × 1° Surface Ocean CO 2 Atlas (SOCAT; Bakker et al., 2016) grid cells contain surface-ocean pCO 2 observations, leading to diverging results in data-sparse regions (Gloege et al., 2021;Ritter et al., 2017).
In contrast to extrapolation-based methods, Ocean Biogeochemistry Models (OBMs; Aumont et al., 2015;Stock et al., 2014) have the ability to resolve the spatiotemporal scales necessary for full attribution of ocean sink variability. OBMs, however, do not typically assimilate physical and biogeochemical observations, which can lead to substantial misfit between the model solution and observations. Until very recently, studies that use data-assimilative OBMs have either been (a) limited to regional studies, for example, the Southern Ocean (Rosso et al., 2017;, or (b) had global-ocean coverage but a short assimilation window (e.g., 2009(e.g., -2011(e.g., in Brix et al., 2015. With the advent of CO 2 -dedicated satellites (e.g., the Observing Carbon Observatory 2 (OCO-2); Crisp et al., 2017) and global carbon-cycle data assimilation systems (e.g., the NASA Carbon Monitoring System Flux (CMS-Flux) project; Ott et al., 2015;Bowman et al., 2017;Liu et al., 2014Liu et al., , 2017Liu et al., , 2021, there is a clear need for global-ocean biogeochemical state estimates that accurately simulate regional patterns in the ocean sink over seasonal to multi-decadal timescales. In fact, OCO-2 observations enabled quantification of spatiotemporal changes in ocean-atmosphere carbon exchange during the 2015-2016 El Niño event (Chatterjee et al., 2017;Patra et al., 2017). Additionally, top-down CMS-Flux estimates have recently been used to inform terrestrial ecosystem dynamics Quetin et al., 2020), which show potential for future applications to ocean biogeochemistry.
Here we use the Estimating the Circulation and Climate of the Ocean-Darwin (ECCO-Darwin) global-ocean biogeochemistry state estimate  to generate a data-constrained DIC budget and investigate how spatiotemporal variability in advection and mixing, air-sea CO 2 flux, and the biological pump have modulated the ocean sink for 1995-2018. This study considers the sum of natural and anthropogenic DIC without, for the moment, separating the two contributions explicitly. ECCO-Darwin assimilates physical (both ocean and sea-ice) and biogeochemical observations, which greatly improves the model's fit to a suite of observations. To our knowledge, this is the first in-depth analysis of the global-ocean DIC budget using an ocean biogeochemistry state estimate. Our results provide novel three-dimensional insight into the bio-chemical-physical decomposition of the DIC pool and trajectory of the ocean sink, priors for global carbon-cycle inversions, and an improved framework for understanding sparse in-situ biogeochemical observations.

ECCO-Darwin: Physical State Estimate
The ECCO-Darwin ocean biogeochemistry state estimate is extensively described in Brix et al. (2015), Manizza et al. (2019), andCarroll et al. (2020). The latest ECCO-Darwin solution (version 5) is based on ocean circulation and physical tracers (temperature, salinity, and sea ice) from the Estimating the Circulation and Climate of the Ocean (ECCO) LLC270 global-ocean and sea-ice data synthesis (Zhang et al., 2018). ECCO LLC270 has nominal 1/3° horizontal grid spacing (∼18 km at high latitudes) and 50 vertical levels. Since the horizontal discretization is insufficient to resolve mesoscale eddies, their impact on large-scale ocean circulation is parameterized using the Redi (1982) and Gent and McWilliams (1990) schemes; vertical mixing is parameterized with the Gaspar et al. (1990) scheme.
Physical observations are assimilated using the adjoint method (i.e., 4-D-Var; Wunsch et al., 2009;Wunsch & Heimbach, 2013), which minimizes a weighted least squares sum of model-data misfit (the cost function) to optimize initial conditions, time-varying surface-ocean boundary conditions, and time-invariant, three-dimensional mixing coefficients for along-isopycnal, cross-isopycnal, and isopycnal thickness diffusivity. Because the initial conditions, surface boundary conditions, and mixing coefficients are estimated as part of the adjoint-method optimization, the ECCO LLC270 ocean circulation estimate has negligible drift and therefore does not require spin-up. The LLC270 circulation estimate is used at each time step (1200 s) to drive an online biogeochemistry and ecosystem model provided by the Massachusetts Institute of Technology Darwin Project Follows et al., 2007) -taken together these components form ECCO-Darwin.

ECCO-Darwin: Biogeochemical State Estimate
The Darwin model includes the cycling of carbon, nitrogen, phosphorus (PO 4 ), iron (Fe), silica (SiO 2 ), oxygen, and alkalinity. Matter is cycled from inorganic nutrients, through living and dead organic matter, and remineralized back to inorganic forms. All particulates are instantly removed once they reach the seafloor as a crude parameterization of sedimentation and to prevent model instability due to the accumulation of particulates at the bottom. Nitrogen cycling explicitly includes nitrate (NO 3 ), nitrite (NO 2 ), and ammonium pools (NH 4 ). Carbonate chemistry is based on the efficient solver of Follows et al. (2006). Air-sea CO 2 flux is computed using the parameterization of Wanninkhof (1992) and forced with apCO 2 from the zonally averaged National Oceanic and Atmospheric Administration Marine Boundary Layer Reference (NOAA MBL) product (Andrews et al., 2014). Atmospheric iron dust deposition at the ocean surface and terrestrial runoff along coastal boundaries is forced using the monthly climatology of Mahowald et al. (2009) and Fekete et al. (2002), respectively. Terrestrial runoff consists of freshwater only and does not include nutrients nor DIC.
The Darwin ecology includes five large-to-small phytoplankton functional types (diatoms, other large eukaryotes, Synechococcus, and low-and highlight adapted Prochlorococcus), along with two zooplankton types that graze preferentially on either large eukaryotes or small picoplankton. Biogeochemical properties (i.e., inorganic carbon and nutrients, phytoplankton and zooplankton, dissolved organic matter, and detrital particles) are treated as prognostic variables and are advected and mixed by the LLC270 ocean circulation.
The biogeochemical model is optimized separately from the circulation using a low-dimensional Green's Functions approach (Menemenlis et al., 2005) to assimilate a variety of biogeochemical observations and to adjust the Darwin initial conditions and six ecological parameters. The mixing coefficients from the adjoint optimization are applied to both the physical and biogeochemical fields. A detailed description of the ECCO-Darwin model setup, observational constraints, and optimization methodology is presented in Carroll et al. (2020).

Model Bias and Drift
An important consideration for the present study is that of model bias and drift. A fully spun-up model will tend to have large bias but negligible drift. Conversely, a model that is initialized close to observations will tend to have small bias but large drift. For the present study, we aim for a compromise between bias and drift that reduces large upper-ocean initialization transients while maintaining the model DIC concentration close to observed profiles in the deep ocean. As shown in Brix et al. (2015), the Green's function approach can simultaneously reduce both bias and drift for biogeochemistry, similar to the adjoint method used for the circulation. Because the Green's Function approach is low dimensional (only a small number of model parameters can be adjusted), some residual model drift is removed by discarding the first 3 years (1992)(1993)(1994) of the ECCO-Darwin simulation in our analysis. Figure 1 shows that the globally integrated ECCO-Darwin air-sea CO 2 flux is consistent with the Jena Carbo-Scope (Rödenbeck et al., 2013) and MPI-SOMFFN (Landschützer et al., 2016) products. The ECCO-Darwin ocean carbon sink trend matches that of the two observation-based estimates, indicating that the ECCO-Darwin global-ocean CO 2 flux has minimal model drift during the 1995-2018 model period. A more detailed region-byregion discussion of ECCO-Darwin CO 2 flux bias and drift is provided in Carroll et al. (2020).

DIC Budget
A key strength of ECCO-Darwin is that all physical and biogeochemical fields are conserved within numerical precision over the entire model period. Therefore, ECCO-Darwin exactly obeys conservation laws and there are no spurious sources and sinks of volume, heat, salt, carbon, etc. This allows for the analysis of a data-constrained, fully closed DIC budget that permits attribution to ocean circulation, air-sea CO 2 flux, and biological processes ( Figure 2).
In the open ocean (i.e., neglecting river and land sources of carbon), the time-evolution of DIC is described by the three-dimensional budget equation, (1) DIC Adv includes the parameterization of transport due to mesoscale eddies from Gent et al. (1995). The diffusion term DIC Dif represents unresolved flows and eddies, for example, convective adjustment, along-isopycnal diffusion from the Redi (1982) scheme, and diapycnal and mixing layer diffusion from the Gaspar et al. (1990) scheme.
Biological processes can be further decomposed as where Org Prod and PIC Prod refer to the biological uptake of DIC for production of organic matter and particulate inorganic matter, respectively. Particulate inorganic matter comes from calcification and is referred to as PIC. PIC Diss is the dissolution of PIC and Org Remin is remineralization of organic matter.
We compute the monthly mean DIC budget in each grid cell during January 1995 to December 2018. The three-dimensional DIC budget is then volume-integrated over the upper 100 m and over the full-depth water column, yielding a mass budget. The 100-m integration depth was chosen because it is similar to the global-mean winter wind-driven mixed layer depth, euphotic zone depth, and is a conventional reference depth for export calculations; potential caveats and sensitivities to this choice will be discussed.
In the ECCO-Darwin model configuration, the injection or removal of surface-ocean freshwater does not impact DIC mass directly, even though it changes surface DIC concentration. The addition or removal of surface freshwater due to precipitation, evaporation, and terrestrial runoff has two consequences. The first is a surface pressure perturbation (or volume change) that is instantaneously distributed globally. The second is a transport of DIC-rich ocean waters away from or toward the location of the surface flux; this effect is captured by the ∇ ⋅ DIC Adv term.
Because exchanges between sea ice and the ocean do not add or remove mass locally, freshwater fluxes due to sea-ice melt and freeze do not transport DIC-rich ocean waters away from or toward the location of the flux. Nevertheless, because the ECCO-Darwin simulation uses the z* rescaled vertical coordinates of Campin et al. (2008), dilution or concentration of DIC due to sea-ice melting or freezing is still, to first order, captured by ∇ ⋅ DIC Adv (Text S1 in Supporting Information S1).  For the remainder of the text, we use the following terminology to represent the DIC budget terms: = DIC tendency, ∇ ⋅ DIC Adv = net advection, ∇ DIC Dif = net diffusion, F CO2 = air-sea CO 2 flux, and DIC Bio = net biology. For the further decomposition of net biology: Org Prod = organic production, PIC Prod = PIC production, PIC Diss = PIC dissolution, and Org Remin = organic remineralization. A summary of all budget terms is shown below in Table 1.

Biome-Scale Analysis
We compute time series of the DIC budget terms in the time-invariant, open-ocean biome regions of . The 5 superbiomes ( Figure 3a) and 17 biomes ( Figure 3b) provide nearly full coverage of the open ocean and are determined from coherence between monthly SST, spring and summer Chl, sea-ice fraction, and climatological maximum MLD for 1998-2010. For simplicity, most of the figures and discussion in this paper pertain to the superbiomes, but we also provide similar analysis for all 17 biomes both in the manuscript and in the Supporting Information for additional detail. These biomes provide a more robust estimate of biogeochemical ocean provinces compared to regions that use arbitrary rectangular boundaries and are a natural choice for clustering carbon cycle dynamics. The biomes are mapped onto the ECCO-Darwin grid using nearest-neighbor interpolation. We note that each biome has a different surface area and hence a different mass-budget volume. Therefore, we normalize by volume  (and compute volume-weighted means) when comparing budget terms across different biomes, resulting in a concentration budget.

Comparison to Observations
We first benchmark ECCO-Darwin climatological DIC (1995DIC ( -2018 with the GLODAPv2 mapped product (Lauvset et al., 2016) and GLODAPv2.2021 in-situ vertical profiles (Lauvset et al., 2021; Figure 4 and Figure  S1 in Supporting Information S1). We note that the previous study of Carroll et al. (2020) extensively evaluated ECCO-Darwin against SOCATv5 surface-ocean pCO 2 observations ; therefore, here we focus only on model-data comparison of DIC. ECCO-Darwin generally reproduces the observed latitudinal gradient in surface-ocean DIC, along with the spatial patterns of low DIC near major river deltas and in the Equatorial, West Pacific, and Indian Oceans ( Figure 4a). In the Greenland, Iceland, and Norwegian Seas, Irminger Sea, Labrador Sea, North Pacific Ocean, and the central Arctic Ocean (where GLODAPv2 lacks observations), ECCO-Darwin generally exhibits higher surface-ocean DIC concentrations compared to the GLODAPv2 mapped product. Examination of the mean (Figure 4c, thin black line) and standard deviation ( Figure 4c, thick blue line) modeldata difference at each depth and in each superbiome reveals a very good fit between model and observations. A maximum standard-deviation difference of 3.9% is reached in the surface level of sNH, resulting from mismatch in locations near sea ice ( Figure 4b, see white outliers in sNH). For all other superbiomes, the maximum standard-deviation difference in the upper 100 m is less than 2.4%. Over the 24-year period of analysis, the majority of DIC variability occurs in the mixed layer and upper ocean, where convection and overturning circulation is vigorous. The standard-deviation model-data mismatch below the 1000-m depth is less than 0.7% in all superbiomes.

Climatological Fields
Having shown that ECCO-Darwin can accurately simulate the surface-ocean patterns and vertical structure of the observed DIC pool, we next examine climatological fields (January 1995 to December 2018) from the DIC mass budget in the upper 100 m ( Figure 5). The DIC tendency (Figure 5a), which is the sum of the budget terms shown in Figures 5b-5e, exhibits a strong latitudinal gradient and shows a general increase of DIC in the tropics, subtropics, and Southern Ocean. The largest gains and losses of DIC occur in the equatorial ocean and along select coastal zones (e.g., the periphery of west Central America and West Africa), respectively. Globally, the largest DIC tendency magnitudes are found in the West Equatorial Pacific Ocean. The Arctic Ocean is a region of weak DIC loss, with negative DIC tendencies concentrated in the central Arctic and offshore of the Kara and Laptev marginal seas.
Net advection results in a gain of DIC in equatorial and coastal upwelling regions and a loss of DIC in the Kuroshio, Gulf Stream, and Agulhas Current systems and their respective extensions (Figure 5b). In upwelling regions, net advection results in a gain of DIC and is dominated by the vertical transport of deep, DIC-rich waters ( Figure 5f). In western boundary currents, net advection results in DIC loss, driven by the poleward transport of tropical, low-DIC waters. The Southern Ocean is also a region of advective DIC gain, albeit weaker than the equatorial and tropical regions, due to the large-scale upwelling (Ekman suction) of DIC-rich waters (Figure 5f). Weaker and more diffuse advective losses of DIC due to downwelling (Ekman pumping) are found in the Northern and Southern Hemisphere subtropical gyres.
We note that the magnitude of the horizontal (not shown) and vertical advection terms ( Figure 5f) are roughly an order of magnitude larger than net advection, demonstrating a substantial level of horizontal-vertical compensation. For example, increased DIC due to horizontal convergence in the Pacific subtropical gyres (not shown) is partially compensated by the divergence of DIC from downwelling ( Figure 5f, blue colors), leading to regions with weak DIC loss. Near the Equator, the vertical upwelling of DIC-rich waters is not fully compensated by horizontal divergence from Ekman transport and poleward flow into the subtropical gyres, leading to a gain of DIC. Globally integrated, net advection provides a DIC gain of 2.9 Pg C year −1 to the upper 100 m.
Net diffusion generally provides a gain of DIC ( Figure 5c) and is dominated by the vertical component due to large vertical DIC gradients in the upper 100 m. Net diffusion hotspots are primarily found in western boundary current systems, where low-DIC waters of tropical origin mix with deeper DIC-rich waters during winter months. At these hotspot locations, positive net diffusion exceeds negative net advection. Therefore, the sum of advection and diffusion, that is, ocean circulation, results in an overall gain of DIC. Globally integrated, net diffusion provides a DIC gain of 3.4 Pg C year −1 . Locally, DIC loss due to diffusion is small and primarily pertains to horizontal diffusion of DIC-rich waters away from upwelling regions. Together, advection and diffusion contribute 6.3 Pg C year −1 to the upper 100 m, which is to a large extent balanced by net biology.
Air-sea CO 2 flux ( Figure 5d) drives CO 2 outgassing (DIC loss) in equatorial regions and offshore of several coastal upwelling zones (e.g., western North America, Central and South America, and northwest Africa). The strongest open-ocean CO 2 outgassing region is located between 2°S and 8°S in the Eastern Equatorial Pacific Ocean. The most vigorous CO 2 uptake (DIC gain) occurs in subpolar and western boundary current regions, where (a) pCO 2 solubility (and hence CO 2 uptake) increases as warm surface waters cool when flowing poleward and (b) strong biological fixation draws down DIC in the euphotic zone. Additionally, poleward-flowing western boundary current waters gain anthropogenic CO 2 as they enter the subpolar gyres, where upwelling exposes waters with low anthropogenic CO 2 concentrations to the surface. In total, air-sea CO 2 flux provides a global-ocean DIC gain of 2.6 Pg C year −1 .
Net biology causes DIC loss in equatorial and high-latitude regions, due to strong organic production that is not fully compensated by organic remineralization in the upper 100 m (Figure 5e and Figure S2 in Supporting Information S1). Near the Equator, increased organic production is driven primarily by the upwelling of nutrients (not shown) and consistent, year-round insolation. In the oligotrophic subtropical gyres, the biological loss of DIC is weak. Globally integrated, net biology results in a DIC loss of 8.6 Pg C year −1 , which has the largest magnitude of all budget terms. Locally, we note that DIC gain occurs where net transport drives convergence of organic material and shallow remineralization outstrips export. The residual between circulation, air-sea CO 2 flux, and net biology (i.e., the DIC tendency) results in a DIC gain of 0.3 Pg C year −1 in the upper 100 m.

Climatological Zonal-Mean Structure
Examination of the zonally averaged mass budget as a function of latitude reveals the intricate pole-to-pole balance between budget terms ( Figure 6). The DIC tendency is roughly an order of magnitude smaller than the other terms, demonstrating substantial compensation in the net budget (Figure 6a, light blue lines). The largest global magnitudes of net advection (dark blue line) and net biology (red line) are found in a narrow equatorial band between 10°S and 5°N, with equatorial upwelling of DIC-rich waters outpacing biological production. Further north, from 5 to 72°N, this pattern reverses and net biology exceeds net advection. Both terms cause a loss of DIC from 30 to 79°N, highlighting the role of strong advective divergence and biological production in the Northern Hemisphere boundary-current systems, Subpolar North Atlantic, and Nordic Seas. In the Southern Hemisphere, from 10 to 70°S, DIC loss from net biology dominates the DIC gain from net advection.
Near the Equator and tropics, net diffusion (light green line) is roughly 8× smaller than net advection. In subtropical and subpolar latitudes (30-60°N and 25-60°S), the magnitude of net diffusion is generally larger than net advection due to mode water formation and intense vertical mixing in boundary currents. These terms tend to oppose each other in the subpolar Northern Hemisphere, with DIC loss from advective divergence being overcompensated by DIC gain from net diffusion. However, in the subpolar Southern Hemisphere, where Southern Ocean upwelling is active, both net advection and net diffusion result in a gain of DIC. Air-sea CO 2 flux (dark green line) generally results in a loss of DIC (CO 2 outgassing) in the equatorial and tropical regions and a gain of DIC (CO 2 uptake) in the subtropics, subpolar regions, and Southern Ocean.
Decomposition of the net biology term (Figure 6b, red line) demonstrates that organic production (light blue line) is partially compensated by organic remineralization (light pink line) at all latitudes in the upper 100 m. The magnitude of both terms is reduced around 30°N and 30°S, where the oligotrophic conditions of the subtropical gyres limit biological productivity and export. PIC production (dark blue line) is roughly two orders of magnitude smaller than organic production, with PIC dissolution (light green line) providing a very small global-ocean gain of DIC.

Seasonal Zonal-Mean Structure
We next examine the seasonal climatology of budget terms as a function of latitude (Figure 7). DIC tendency exhibits clear seasonality in both hemispheres, with a stronger seasonal cycle in the Northern Hemisphere. Throughout the equatorial band, intricate patterns in net advection, which shift latitudinally throughout the year, drive the largest advective seasonal cycle in the global ocean. In the Arctic Ocean and periphery of Antarctica, the formation and melt of sea ice causes advective divergence and convergence of DIC, respectively; this effect is strongest in the Arctic Ocean during boreal summer. Seasonality in net diffusion results from increased vertical mixing during winter and spring. This occurs primarily in boundary current regions and throughout mode water formation regions in the Northern North Atlantic and Southern Ocean (Figure 7b). This is due to intense vertical mixing in the Gulf Stream, Kuroshio, Agulhas, and East Australian Current systems.
Both air-sea CO 2 flux and net biology have a seasonal structure that strongly varies with latitude. For air-sea CO 2 flux, this latitudinal shift is driven by a combination of physical and biogeochemical processes that cause seasonal variations in surface-ocean pCO 2 solubility, wind speed, and mixed layer depth. Air-sea CO 2 flux seasonally alternates between CO 2 outgassing (DIC loss from summer-fall) and uptake (DIC gain from winter-spring) in both hemispheres, with year-round outgassing near the Equator (Figure 7b). Net biology produces a consistent loss of DIC at the Equator, with spring blooms resulting in seasonal DIC loss in both hemispheres.

Seasonal Biome-Scale Analysis
We now zoom into the biome scale to examine the amplitude of the climatological seasonal cycle for each budget term ( Figure 8 and Table S1 in Supporting Information S1). Here the mass budget is normalized by volume to yield concentration, allowing for a comparison between biomes that is not biased by surface area. For all superbiomes, sNH has the largest DIC tendency amplitude (∼0.2 mol m −3 year −1 ), followed by sSH (Figure 8a). Seasonality in sNH results primarily from (in decreasing order) net biology, net advection, and air-sea CO 2 flux (Figure 8b). Examining the individual biomes that comprise sNH, we find that seasonality in high-latitude NP1 and NA1 is largely modulated by net advection, with net biology playing a smaller role (see Figure S3 in Supporting Information S1). Here, summertime sea-ice melt dilutes DIC-rich waters, leading to divergence in net advection. Moving to the subpolar, seasonally stratified NP2 and NA2, the seasonal cycle is principally driven by strong net biology and air-sea CO 2 flux. Further south, in the subtropical seasonally stratified NP3 and NA3 biomes, which include the Gulf Stream and Kuroshio Current, net biology is weaker and net diffusion increases due to elevated vertical mixing.
In sSH, which is the Southern Ocean superbiome, DIC tendency is dominated by net biology (SO1 and SO2) and net advection along the periphery of Antarctica (SO3). Here, SO3 is strongly impacted by sea-ice melt during austral summer ( Figure S3 in Supporting Information S1). In equatorial sEQ, seasonality in DIC tendency is driven by time variability in net advection (WPE, EPE, and AE), which reaches local minima during the equinoxes; net biology and air-sea CO 2 flux have a weak seasonal cycle, with net diffusion having negligible seasonality ( Figure 8b). For sSS, the seasonal cycle results from a combination of spring blooms (SP and SA), winter CO 2 outgassing and summer uptake (SP, SA, and IO), and net advection in the Indian Ocean (IO). sNS, which contains permanently stratified NP4 and NA4, has the smallest DIC tendency amplitude of all superbiomes.  Table S1 in Supporting Information S1. (b) Seasonal climatology of concentration budget terms in the superbiomes; dashed black line shows the zero value. See Figure S3 in Supporting Information S1 for all individual biomes.

Global-Ocean Upper-100-m Time Series
Having shown climatological and seasonal results, we next examine global-ocean interannual variability (IAV) in the upper 100 m and full-depth DIC mass pool ( Figure 9). DIC content in the upper 100 m exhibits strong IAV (Figure 9a), which is associated with the ENSO signal. ENSO events cause a sharp loss and subsequent gain in global-ocean DIC mass during alternating El Niño and La Niña years, respectively (e.g. , 1997-1998 and 2015-2016). Over the 24-year model period, the 1997-1998 ENSO caused the largest global year-to-year change in upper-ocean DIC mass, varying by 2.1 Pg C over two years. A least squares fit of annual-mean DIC during 1995-2018 highlights the linear increasing trend in the upper ocean (Figure 9a, dashed black line). Over annual timescales, the upper-100-m DIC pool tracks the atmospheric CO 2 growth rate, which to first-order is linear over our model simulation. During 1995-2018, the upper-100-m DIC pool increased by 8.1 Pg C. Examination of the associated mass budget reveals that IAV in the upper-100-m DIC pool is primarily driven by net advection and is due to ENSO (Figure 9b). Air-sea CO 2 flux and net diffusion have opposing trends, as the secular increase in ocean CO 2 uptake reduces upper-ocean carbon gradients and weakens vertical mixing of DIC. Over the 24-year model period, biological DIC loss in the upper 100m weakened by roughly 6.7%.

Global-Ocean Full-Depth Time Series
For the full-depth DIC pool (Figure 9c), IAV from ENSO events has a negligible impact on the large ocean reservoir, with a small seasonal cycle driving subannual variability (gray line). Because of increasing atmospheric CO 2 concentration, the global-ocean DIC pool cumulatively sequesters increasing atmospheric CO 2 , resulting in a slightly nonlinear increase of DIC mass. Over the model period, the global-ocean DIC pool increased from 37,455 to 37,519 Pg C (or Gt C), a difference of roughly 64 Pg C (∼2.7 Pg C year −1 ). In the corresponding fulldepth budget (Figure 9d), the secular increase in DIC tendency, which increases by roughly 64% over the model period, is driven primarily by the linear increase in air-sea CO 2 flux. Full-depth net biology provides a small contribution that is generally centered around zero (time mean of 0.1 Pg C year −1 ).
We note that net advection and diffusion fully compensate when integrated across the full-depth global ocean, that is, circulation only redistributes DIC and does not result in a net gain or loss. Also note that PIC/POC removed at the model seafloor is added to the full-depth DIC budget in Figures 9c and 9d. This is because in the real ocean, particulate carbon is expected to either be remineralized and dissolved before being buried or to be replenished by riverine fluxes, processes that are not currently represented in ECCO-Darwin.

Biome-Scale Time Series
Having shown the global-ocean results, we next focus again on the biome scale to examine regional IAV. Annual-mean time series of concentration budget terms (Figure 10), taken in the upper 100 m of all biomes, demonstrates that DIC tendency (light blue line) closely tracks IAV in net advection (dark blue line). Over the 24-year model period, DIC tendency is generally positive (DIC gain) when averaged across all superbiomes, only becoming slightly negative (DIC loss) during major ENSO events. sNH exhibits negative net advection during the entire model period; the magnitude of this DIC loss results from consistent advective divergence in NA2 and NA3, with associated IAV driven by ice-influenced NP1 and NA1 ( Figure S4 in Supporting Information S1). In sNS, DIC tendency is also dominated by sharp year-to-year variability in net advection, with IAV in air-sea CO 2 flux and net biology playing a secondary role.
Equatorial sEQ shows the clear influence of ENSO events and other long-term fluctuations. Here net advection undergoes a sharp reduction during El Niño years due to arrested equatorial upwelling in WPE and EPE ( Figure  S4 in Supporting Information S1). The reduction in (a) net advection of DIC-rich waters to the mixed layer and (b) nutrients to the euphotic zone during El Niño (not shown) causes a corresponding decrease in CO 2 outgassing and biological uptake, respectively. We note that the very strong decrease in net advection causes only a minor reduction in net outgassing, due to compensation from thermal effects (i.e., warmer waters cause more outgassing of natural CO 2 ) and anthropogenic CO 2 uptake (reduced upwelling results in less anthropogenic CO 2 uptake). In WPE, the 2015 El Niño reduced CO 2 outgassing enough that short-term uptake was achieved. Interestingly,  Figure S4 in Supporting Information S1 for all individual biomes.
net diffusion in sEQ does not exhibit the strong secular decline seen in the other superbiomes; this is because CO 2 outgassing helps maintain the vertical structure of the carbocline, holding net diffusion steady. We find that DIC tendency in sSS is also impacted by major ENSO events (e.g., 1997-1998 and 2015-2016), highlighting the reduced upwelling and poleward advection occurring in both sEQ and sSS.
In Southern Ocean sSH, IAV is driven by a combination of net advection, air-sea CO 2 flux, and net biology.
Here IAV predominately occurs in SO2 and SO3 ( Figure S4 in Supporting Information S1), which are located closer to the Antarctic continent. The sum of terms with DIC gain generally exceeds the loss from net biology, except during 1996, 2002, 2009, and 2016; years when net advection in the ice-dominated SO3 decreased sharply ( Figure S4 in Supporting Information S1), leading to a negative DIC tendency.

Biome-Scale Magnitudes and Trends
Finally, we evaluate biome-scale budget magnitudes and trends, computed from annual-mean concentration budget time series in the upper 100 m (Figure 11 and Table S2 in Supporting Information S1). The small timemean DIC tendency values, which are roughly an order of magnitude smaller than the other terms, again highlight the substantial compensation between bio-chemical-physical processes (Figure 11a).
For all superbiomes, equatorial sEQ has the strongest time-mean net advection, resulting from vigorous upwelling of DIC-rich waters along the Equator in EPE and AE. Maximum net advection (DIC gain of ∼0.05 mol m −3 year −1 ) is obtained in EPE, which is located in the eastern Equatorial Pacific Ocean. In contrast, Northern Hemisphere subpolar and subtropical seasonally stratified biomes (NP2, NP3, NA2, and NA3) have negative time-mean net advection, primarily driven by the advective divergence of DIC in boundary current systems and their poleward-flowing extensions. Ice-influenced biomes NP1 and NA1 also exhibit negative net advection, caused by the dilution of DIC-rich waters from sea-ice melt during summer. All superbiomes have positive timemean net diffusion, with sNH having the largest diffusive gain of DIC (∼0.01 mol m −3 year −1 ), due to elevated vertical mixing in the subpolar and subtropical seasonally stratified biomes (NP2, NP3, NA2, and NA3). Positive  Table S2 in Supporting Information S1.
air-sea CO 2 flux (uptake) is found in all superbiomes except for equatorial sEQ, which has a negative time-mean flux (outgassing). A maximum superbiome air-sea CO 2 flux of ∼0.02 mol m −3 year −1 occurs in sNH, driven by strong CO 2 uptake in NP3 and NA2. For all individual biomes, maximum and minimum air-sea CO 2 fluxes occur in North Atlantic NA2 and equatorial EPE, respectively. Time-mean net biology is negative (DIC loss) in all superbiomes, with the largest magnitudes (in decreasing order) occurring in sEQ, sNH, and sSH. For all individual biomes, the strongest biological uptake occurs in Southern Ocean SO1 (∼0.04 mol m −3 year −1 ).
The majority of statistically significant trends (95% confidence interval) occur in net diffusion, air-sea CO 2 flux, and net biology terms (Figure 11b, white crosses show p-values ≤ 0.05). For all biomes, significant trends in DIC tendency occur only in South Atlantic SA and Southern Ocean SO1, with both biomes having a positive trend. Significant trends in net advection are obtained in 7/17 biomes, with Indian Ocean IO and Southern Ocean SO3 having the only negative trends. All superbiomes and biomes have statistically significant trends in net diffusion, with negative trends in all biomes except for the ice-influenced NA1. The majority of biomes (15/17) demonstrate significant trends in air-sea CO 2 flux, with Northern Hemisphere ice-influenced NP1 and NA1 having the only significant negative trends. Significant trends in net biology occur in sNH and sSH, and in roughly 70% of biomes (12/17). Only NP4 and SP have significant negative trends in net biology (i.e., increasing biological DIC loss).

Overview
Data-based reconstructions provide increasingly accurate estimates of the ocean sink for atmospheric CO 2 . These reconstructions, however, do not fully explain the mechanisms that drive space-time variability in the sink. Alternatively, hindcast OBMs provide a tool for exploring the individual processes that control ocean carbon variability; however, most of these models are not constrained by observations. Using the ECCO-Darwin global-ocean biogeochemistry state estimate, we have combined the above two methods and produced a multi-decadal , data-constrained DIC budget that fully attributes ocean sink variability to three-dimensional physical, chemical, and biological drivers.
The fully closed decomposition of DIC budget terms is a powerful tool for explaining observed DIC variability. In the upper ocean (top 100 m in this study), advective and diffusive processes (i.e., circulation) provide the largest gain of DIC, with net biology resulting in the largest loss. When considering the volume-integrated global ocean, air-sea CO 2 flux provides the largest DIC gain (time mean and standard deviation of ∼2.6 ± 0.5 Pg C year −1 , respectively). Globally, interannual variability is dominated by vertical advection in equatorial regions, with ENSO causing the largest year-to-year change in upper-ocean DIC mass. In terms of the altered modern vertical DIC gradient, we note that a natural DIC flux, which has similar magnitude to the loss by net biology, is brought to the surface by advective and diffusive processes (i.e., ocean circulation). However, there is also a residual downward anthropogenic DIC flux from air-sea CO 2 exchange. Assuming the imbalance between circulation and net biology roughly reflects the net downward anthropogenic DIC flux (2.3 Pg C year −1 ), there would be a physical anthropogenic carbon removal of 55 Pg C to depths greater than 100 m over the 24-year model period. This accounts for 85% of the full-depth DIC mass increase over the same period, with ∼8 Pg C remaining in the upper 100 m. We next summarize the time-mean, blended patterns of DIC gains and losses across the global ocean ( Figure 12).

Circulation
In the upper 100-m, ocean circulation (Figure 12a, purple colors), which is the sum of net advection and net diffusion, is the primary gain of DIC in equatorial regions, western boundary currents, and along eastern boundary current upwelling zones (e.g., California, Peru, Benguela, and Canary Current systems). While we note that diffusion represents unresolved turbulent processes and the split between advection and diffusion depends on model configuration and the choice of the 100 m vertical integration depth, we hypothesize that the combined circulation component is robust. For circulation in western boundary currents, diffusion exceeds advection, while in the equatorial and upwelling regions advection dominates circulation (Figures 5b and 5c).
When examining the full-depth ocean (Figure 12b), gains of DIC from intense swaths of circulation are reduced in size, due to substantial compensation of advective and diffusive processes throughout the water column. In the upper 100 m, circulation results in a weak loss of DIC only in few select regions (e.g., subpolar North Atlantic and Indian Oceans). For the full-depth ocean, DIC loss from circulation is substantial outside of equatorial regions, highlighting the large-scale extratropical divergence required to balance equatorial convergence of DIC-rich waters.

Air-sea CO 2 Flux
Ocean-atmosphere exchange of natural and anthropogenic CO 2 (Figure 12, orange colors), which is the same in the upper 100 m and full-depth ocean, results in a gain of DIC in western boundary currents, subpolar regions, select subtropical regions (e.g., Pacific and Indian Oceans), and along the Subantarctic Front of the Antarctic Circumpolar Current (ACC). For subtropical western boundary currents, warm waters are transported poleward, which leads to dramatic cooling and a subsequent increase in CO 2 solubility. Meanwhile, subpolar boundary currents bring cold waters equatorward, which then warm and hence outgas CO 2 . Additionally, subtropical western boundary currents are regions exhibiting strong biological loss of DIC, which allows for larger air-sea pCO 2 gradients and hence increased CO 2 uptake. Air-sea CO 2 flux drives a loss of DIC in equatorial regions, where upwelling of DIC-rich waters from circulation results in intense CO 2 outgassing, and in broad swaths offshore of several eastern boundary currents (e.g., California and Peru Current systems).

Biological Processes
In the upper 100 m, biological processes (Figure 12a, green colors) result in a loss of DIC due to fixation of carbon in the euphotic zone and export of particulate carbon to depth. Hotspots of DIC loss from upper-ocean biology directly overlap with circulation-driven DIC gain as upwelling and vertical mixing of carbon also delivers vital nutrients to the euphotic zone, fueling biological productivity. The upwelling that causes strong productivity also results in horizontal divergence away from highly productive regions. In the full-depth water column, subtropical regions are locations of DIC gain (Figure 12b). This is primarily due to lateral convergence of dissolved organic carbon (DOC), which is then remineralized away from regions where it was produced, that is, regions of high biological productivity. While particulate organic carbon (POC) sinks and is remineralized below the export region, DOC remains in suspension and is advected laterally to the center of the subtropical gyres. Compared to the upper 100 m, full-depth biological loss is reduced in both intensity and surface area due to DOC being transported and remineralized away from regions of high productivity.

Global-Ocean DIC Mass
The simulated global-ocean DIC mass from ECCO-Darwin (time mean of 37,484 Pg C) agrees well with previous broad-scale estimates of the inorganic carbon pool (37,400 Pg C in Falkowski et al., 2000) and recent databased reconstructions (37,214 ± 195 Pg C in Keppler et al., 2020). Additionally, our results suggest a full-depth DIC gain of ∼64 Pg C between 1995 and 2018, with air-sea CO 2 flux contributing a total of 63 Pg C and the biological pump contributing 1.4 Pg C (Figure 9d). In the upper 100-m, we estimate a smaller ocean sink increase of 8.1 Pg C over the model period, accounting for roughly 13% of the full-depth increase. In the same depth range, net biology result in a net community production (NCP; Sarmiento & Gruber, 2006) of 8.6 Pg C year −1 , which is in line with previous observation-based (Lee, 2001;Li & Cassar, 2016) and data-assimilative model (9 Pg C year −1 ; DeVries & Weber, 2017) estimates, lending confidence to our estimates of global-ocean productivity. Furthermore, our ocean CO 2 sink trajectory is generally consistent with reconstructions derived from observations (Gruber et al., 2019;Khatiwala et al., 2009) These results highlight the utility of using a global-ocean biogeochemistry state estimate, such as ECCO-Darwin, to: (a) complement data-based reconstructions and quantify the individual mechanisms that modulate ocean sink variability (McKinley et al., 2017); and (b) provide critical constraints on the long-term impact of the terrestrial biosphere on atmospheric CO 2 concentrations (Friedlingstein et al., 2019). Additionally, our results from the 1995-2007 period suggest that 53% of the atmospheric CO 2 uptake between 1995 and 2018 occurred after 2007, with the later 11 years having a time-mean CO 2 exchange rate of roughly 3 Pg C year −1 . Finally, we note that the Green's Functions approach described in Carroll et al. (2020), combined with sensitivity experiments that separate the global-ocean DIC inventory into natural and anthropogenic components (DeVries, 2014;Sabine et al., 2004), could be used to gain a more complete understanding of natural and forced variability in the ocean sink (McKinley et al., 2020).

Interannual Variability and ENSO
Globally, we find that the long-term trend and IAV in upper-ocean DIC is driven principally by (a) the growth rate of atmospheric CO 2 and (b) major ENSO events, respectively ( Figure 9a); results that are supported by previous model-based studies (McKinley et al., 2004studies (McKinley et al., , 2020Resplandy et al., 2015). We also note that our climatological DIC tendency of 0.3 Pg C year −1 generally matches observation-based estimates of ∼1 μmol kg −1 year −1 (Gregor & Gruber, 2021). In the Pacific Ocean between 15°S and 8°N, vigorous outgassing of CO 2 drives a consistent loss of DIC during ENSO-neutral years ( Figure 6 and Figure S4 in Supporting Information S1). This CO 2 outgassing is stronger on the eastern side of the basin (Figure 5d), where ocean-atmosphere pCO 2 gradients are large Landschützer et al., 2016;Sutton et al., 2014). However, during major El Niño events, equatorial upwelling of DIC-rich waters is substantially arrested and net advection decreases, causing an associated decline in CO 2 outgassing and biological processes ( Figure 10, see sEQ). The changes in carbon exchange were large enough during the 2015-2016 El Niño event that they were observed by OCO-2 (Chatterjee et al., 2017;Patra et al., 2017). During La Niña, equatorial upwelling is reinvigorated, with CO 2 outgassing and biological productivity rebounding sharply.
Examination of the Pacific Ocean DIC budget in Hovmöller space (Figure 13), demonstrates the wide range of bio-physical ENSO diversity (Capotondi et al., 2015;Gierach et al., 2012) captured by our simulation. We find that during the major 1997-1998 and 2015-2016 El Niño events (vertical gray lines), net advection along the Equator changes sign and transitions from a DIC gain to a weak loss. Time-depth Hovmöller diagrams, taken in WPE and EPE, highlight the asymmetrical west-to-east response of the equatorial carbocline and budget terms to diverse ENSO events ( Figure S5 in Supporting Information S1).
In west equatorial WPE, the carbocline is deeper and is capped by low-DIC waters associated with the warm pool (Feely et al., 2002). During the large 1997-1998 and 2015-2016 ENSO events, WPE exhibits increased upwelling of DIC-rich waters and biological production; this reduces outgassing to near zero ( Figure S5a in Supporting Information S1). An opposing set of changes occurs in east equatorial EPE during ENSO, with El Niño conditions resulting in a sharp loss of DIC in the upper 100 m, which corresponds to arrested equatorial upwelling, reduced outgassing of CO 2 , and a decrease in biological production and deep remineralization ( Figure  S5b in Supporting Information S1). Additionally, our results show that advective ENSO signals along the Equator are propagated poleward by Ekman transport and subsequent flow into the subtropical gyres ( Figure 10, see net advection in sEQ and sSS).

Observed Trend or Model Drift?
There are two standard ways of initializing an ocean biogeochemistry model. The first is to start the simulation close to observations, for example, a climatology, which can result in a large model drift over time. The second is to start the simulation after a long spin-up, which can result in a large model bias relative to observations. The third approach, and the one that we adopt here for both the physics and biogeochemistry, is to use data assimilation to obtain a solution that simultaneously minimizes model bias and drift, according to an objective cost function. Rather than spin-up the model for thousands of years to a preindustrial steady-state before applying an anthropogenic CO 2 perturbation, we aim to obtain 1995 initial conditions that have small bias compared to observations at depth (see Figure 4c and Figure S1 in Supporting Information S1) while at the same time minimizing the upper-ocean initialization shock. Because of the low-dimensionality of the biogeochemical estimation (i.e., small number of control parameters), we further reduce initialization shock by discarding the first 3 years (1992)(1993)(1994) of the simulation ( Figure S6 in Supporting Information S1). The dashed lines in Figure 9 show simulated trends during 1995-2018. A key question is whether these trends are real or result from residual model drift. Although we have used data assimilation to minimize model bias and drift, we acknowledge that the drift of some of the budget terms, for example, the net diffusion or net biology terms in Figures 9b and 10, may be a residual model spin-up issue as opposed to a real signal.
The net diffusion term for the upper 100 m of the global ocean (Figure 9b, light green line) has a 1995-2018 trend of −72 TgC year −2 . This trend is primarily caused by the increasing anthropogenic perturbation of atmospheric CO 2 , which tends to decrease the vertical gradient of total DIC in the upper ocean. We have verified this via a sensitivity experiment that reduced apCO 2 , maintaining the observed 1995 seasonal cycle for the complete 1995-2018 period. In this sensitivity experiment with lower apCO 2 , the diffusion term contribution to upper-100-m DIC mass is higher, that is, it decreases more slowly during the 1995-2018 integration period, as would be the case when a sharper vertical DIC gradient is maintained. Note that the decreasing net diffusion trend is accompanied by increasing CO 2 flux trend (Figures 9b and 10, dark green line). Also note that the Equatorial region, where there is little anthropogenic CO 2 uptake, has no trend in the diffusion term of the upper 100 m (see sEQ in Figure Figure S7 in Supporting Information S1) indicates that model drift may, in part, be the cause for this trend. However, the biology trend is small relative to interannual variability and it is only noticeable for globally integrated quantities. In all cases, there are insufficient observations to verify globally integrated, model-predicted trends. These model predictions will need to be confirmed or invalidated by future observational or modeling studies.
We also stress that the choice of vertical integration depth for the DIC mass budget affects compensation and hence the relative importance of budget terms. The Southern Ocean DIC budget study of Rosso et al. (2017) used a depth of 650 m, which was chosen to span the deepest modeled mixed layer. We find that using a shallower integration depth of 100 m results in increased net diffusion (due to atmospheric CO 2 uptake), a term which was shown to have a small impact in Rosso et al. (2017). Globally, we find that depth-integrated net advection exceeds net diffusion at depths greater than 127 m (not shown). The points raised above should be considered when interpreting our upper-ocean budget results.

Model Caveats and Potential Improvements
In addition to being a powerful tool for explaining observed DIC variability, the decomposition of DIC budget terms can also be used to identify model deficiencies and their causes. A partial list of known model deficiencies are enumerated below as a focus for future work. Our treatment of freshwater runoff along the coastal zone, and from the Greenland and Antarctic Ice Sheets, is coarsely represented and does not capture IAV (Fekete et al., 2002). Furthermore, we do not include biogeochemical runoff (i.e., DIC and nutrients, Mayorga et al., 2010), which may be critical for accurate simulation of Arctic Ocean biogeochemistry (Le Fouest et al., 2013Fouest et al., , 2015. Future efforts should incorporate time-varying, point-source global freshwater (Feng et al., 2021) and biogeochemical runoff (Resplandy et al., 2018) to improve the simulated ocean sink in coastal (Chen et al., 2013) and high-latitude regions (Hopwood et al., 2018(Hopwood et al., , 2020Wadham et al., 2019). Finally, ECCO-Darwin does not contain an explicit representation of bottom sediment dynamics. To balance the input of riverine carbon and to help address large uncertainties in carbon burial rates (Dunne et al., 2007), future efforts should include a realistic representation of burial and sediment-water fluxes, that is, add a bottom sediment model (Sulpis et al., 2021) and optimize dissolution, remineralization, and particle-sinking rates.
The ECCO-Darwin biogeochemical optimization also relies on a low-dimensional (i.e., small number of control variables) Green's Functions approach to minimize model-data misfit Menemenlis et al., 2005). The incorporation of an adjoint-based optimization method, as done in Verdy and Mazloff (2017) and Rosso et al. (2017) for the Southern Ocean, would allow for a larger number of control variables and joint physical-biogeochemical optimization, which could further improve the fit of simulated DIC to observations. Future work will incorporate the Darwin spectral light and radiative transfer package (Dutkiewicz et al., , 2018(Dutkiewicz et al., , 2019, which will allow for the assimilation of satellite ocean color data. Combined with in-situ particle flux observations (Estapa et al., 2019;Siegel et al., 2016), this would: (a) provide critical constraints on simulated phytoplankton dynamics, export, and remineralization rates; and (b) improve the model representation of climate-driven variability in biological export, for example, decline in primary productivity due to increased stratification (Lozier et al., 2011).
Finally, we stress that further model-data uncertainty analysis is required to confirm the multi-decadal trends evident in our budget analysis (Figure 10b). We note, however, that our use of data-constrained circulation estimates from ECCO lends credibility to the simulated space-time variability in advective and diffusive budget terms. Furthermore, ECCO-Darwin surface-ocean pCO 2 and air-sea CO 2 flux has been extensively evaluated against leading data-based reconstructions (Landschützer et al., 2016;Rödenbeck et al., 2015) and shows broadscale consistency across most biomes, particularly in the subtropical and equatorial regions . As mentioned previously, further efforts and comparison with observations are needed to confirm that trends in our simulated biology are driven by natural variability. Additionally, we acknowledge that the nominal 1/3° horizontal grid spacing of ECCO-Darwin does not permit explicit representation of mesoscale eddies and submesoscale processes. As it becomes computationally practical, we expect that increased model resolution will allow for representation of fine-scale, bio-physical interactions (Mahadevan, 2016;Whitt et al., 2019).

Concluding Remarks
We have produced a data-constrained, global-ocean DIC budget that permits attribution of space-time variability in the ocean sink to three-dimensional ocean circulation, ocean-atmosphere CO 2 exchange, and biological processes. The data-constrained DIC budget presented here represents a significant advance and clearly demonstrates the intricate spatial diversity and compensation in ocean sink mechanisms. These results provide much-needed context for sparse in-situ biogeochemical observations and can be used to inform data-based reconstructions of the global-ocean DIC pool and ocean sink trajectory. The ECCO-Darwin solution, the budget terms, and all software used to generate the simulation and diagnostics have been made available to the science community and are expected to enable further global and regional biogeochemical studies. As (a) the fidelity of ECCO physical-ocean circulation estimates continue to improve and extend in time and (b) upper-ocean biogeochemical observations become more numerous (e.g., from BGC-Argo floats), we anticipate that our data-constrained DIC budget will be an ever more useful tool for ocean carbon sink, ecosystem, and climate-related studies. Ultimately, this work supports the future development of a fully coupled ocean-atmosphere carbon state estimate.

Data Availability Statement
ECCO-Darwin model output and DIC budget diagnostics are available at the ECCO Data Portal: http://data.nas. nasa.gov/ecco/. Model code and platform-independent instructions for running ECCO-Darwin simulations are available at: https://zenodo.org/record/6091603. A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004), with grants from ocean biology and biogeochemistry, physical oceanography, modeling, analysis, and prediction, and interdisciplinary studies programs. High-end computing resources were provided by the NASA Advanced Supercomputing (NAS) Division of the Ames Research Center. Government sponsorship acknowledged. MM acknowledges US National Science Foundation (NSF) awards PLR-1425989 and OPP-1936222. JML acknowledges NSF award PIRE-1545859 and the Simons Collaboration on Computational Biogeochemical Modeling of Marine Ecosystems (CBIOMES) award 549931. We thank Tim DeVries and an anonymous reviewer for their thoughtful and constructive criticism of this work.